%% Colebrook2
function [Pz]=Colebrook2(Tpj,Ppj,d,Z_s,relative_den)
%% lcgc
load lcgc0.mat
StationPoint=[9 24];
lcgc0(:,2)=lcgc0(:,2)-lcgc0(1,2);
F_area1=zeros(StationPoint(1)-1,1);
for lc=1:StationPoint(1)-1
    F_area1(lc)=(lcgc0(lc,2)+lcgc0(lc+1,2))*(lcgc0(lc+1,1)-lcgc0(lc,1));
end
s_z=[lcgc0(StationPoint(1),2) lcgc0(StationPoint(2),2)-lcgc0(StationPoint(1),2)];

lcgc0(:,2)=lcgc0(:,2)-lcgc0(StationPoint(1),2);
F_area2=zeros(StationPoint(2)-StationPoint(1),1);
for lc=StationPoint(1):StationPoint(2)-1
    F_area2(lc-StationPoint(1)+1)=(lcgc0(lc,2)+lcgc0(lc+1,2))*(lcgc0(lc+1,1)-lcgc0(lc,1));
end
F=[sum(F_area1) sum(F_area2)];

%% data
L=[60 100];
Q=[6 1.7143] ; % 10^6Nm3/d
C0=104.9876;

Z=mean(Z_s,2);
d_c=zeros(3,2);
d_c=[d(:,2) d(:,5)];

Tpj=Tpj+273.15;
d_c=d_c/1000;

Pz=zeros(size(d_c));
P1=zeros(size(d_c));
P1(:,1)=[5 5 5]';
P1(:,1)=P1(:,1)+101.325*10^-3;
%% loops
for Station =1:2
    for Plan=1:3
        [lambda,a_coeff]=lambda_cal(Tpj(Plan)-273.15,Ppj(Plan),d_c(Plan,Station)*1000,Q(Station)*10^6/(24*3600),relative_den);
        Pz(Plan,Station)=((P1(Plan,Station)^2-...
            Q(Station)^2/(C0^2*d_c(Plan,Station)^5)*lambda*Z(Plan)*Tpj(Plan)*relative_den*L(Station)*(1+a_coeff*F(Station)/(2*L(Station))))...
            /(1+a_coeff*s_z(Station)))^0.5;
        
        P1(Plan,2)=Pz(Plan,1);
    end
end

end
